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Abstract 

We propose a simple random process inducing various types of random graphs and the scale free random 
graphs among others. The model is of a threshold nature and differs from the preferential attachment 
approach discussed in the literature before. 

The degree statistics of a random graph in our model is governed by the control parameter rj stirring 
the pure exponential statistics for the degree distribution (at rj — 0, when a threshold is changed each time 
a new edge added to the graph) to a power law (at r\ = 1, when the threshold is frozen). The exponent 7 
characterizing the power law can vary in the wide range 7 £ (1, 00) and can be tuned in different values 7i n 
and 7out for in-degrees and out-degrees probability distributions independently. For the intermediate values 
of rj, the decay rate is mixed. 

Taking different statistics for the threshold changes, one obtains dissimilar asymptotic profiles for the 
degree distribution having, in general, nothing to do with power laws at 77 = 1, but still uniformly exponential 
at rj = 0. 
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1 Introduction 



A wide variety of systems in biology, communication technology, sociology and economics are 
best described as complex networks of agents linked by various physical and informational links. 
Clearly, the architecture of such networks is crucially important for the vitality and performance 
of these systems. It is generally recognized that a thorough investigation of possible mechanisms 
which determine the topology of complex "real- world" networks is necessary for understanding 
their behavior. Numerous studies of the Internet and World Wide Web (WWW) ||-[|], 

the web of human sexual contacts || and movie actors collaboration networks ]7|-[[L1[|, national 



power-grids |H] and phone-call networks [T2 

ni 



protein folding 



and cellular networks 



2^|, etc. have revealed, among other facts, that the probability that a node of these 



networks has k connections follows a power law 

P(k) oc fc" 



(1) 



over a large range of k, with an exponent 7 that ranges between 1 and 3 depending on the system 
|25fl . Such networks are called scale- free @. Many other "real- world" networks display an 
inherently uniform exponential topology that still deviates from a Poisson distribution expected 
for usual random graphs G(N,p). 

Recently, a model generating a scale-free network has been suggested in [^J. It establishes 
a random graph process in which vertices are added to the graph one at a time and connected 
to a fixed number of earlier vertices, selected with probabilities proportional to their degrees. 
This preferential attachment assumption is based on the idea that a new site is more likely to 
link to existing sites which are "popular" at time the site is added. Numerical simulations for 
the above process reported in [|| [26|] show that after many iterations the proportion of vertices 
with degree k indeed obeys the power law (|l|) with the exponent 7 = 3. In []27|, the power law 
statistics for P(k) in this model has been proven analytically, and it has been justified that 
7 = 3. 

Probably, a complex network formed by the citation patterns of scientific publications where 
the nodes standing for published articles and a directed edge representing a reference to a 
previously published article exhibiting the power law with the exponent 7 = 3 [ 2£| , can be well 
fitted by the above process as discussed in P, f 



27 . However, in other "real- world" networks 



having degree sequences matching a power law these exponents usually deviate from 3. For 
example, 7 = 2.3 ±0.1 for the movie actors collaboration network J7|, RL 7 = 2.1 for a graph 



of long distance telephone calls made during a single day ]T2], [U| , etc. To fit these deviations, 
one assumes that the preferential attachment for a real network is nonlinear [TJL pl|, i.e. the 
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probability that a node i attaches to node j can be characterized by some function IT(fc) also 
following a power law with some index a. 

The effect of nonlinear preferential attachment H(k) on the network dynamics and topology 
was studied in [|l(|. They demonstrated that the scale free network topology holds only if 
H(k) is asymptotically linear and is destroyed otherwise. This way the exponent of the degree 
distribution P(k) can be tuned to any value between 2 and oo [CT. 

Nevertheless, the detailed studies of various "real- world" networks still challenge to the 
preferential attachment approach |J. Some ecological networks such as food webs quantifying 
the interactions between species []l7|], some networks in linguistics fl^l , a network representing 
an aggregation of the WWW domains [18], etc. exhibit a scale free topology characterized by 
the indices 1 < 7 < 2. Moreover, for some networks, the in-degree and out-degree probability 
distributions are dissimilar: the WWW-topology where the nodes are the webpages and the 
edges are the hyperlinks that point from one document to another meets the power laws with 
7i n = 2.1 and 7 out = 2.72 [29], the distributions of the total number of sexual partners in a 
single year indicate the power law decay with 7 = 2.54 ± 0.2 for females and 7 = 2.31 ± 0.2 for 
males ||. A recent study |3(| extends the results on the citation patterns (incoming links) of 
scientific publications (7^ = 3, p8j) to the outgoing degree distribution obtaining that it has 
an exponential tail (not a Poisson distribution). This challenge stimulates a search for the new 
models inducing graphs of the scale free topology. 

A random procedure which we use to generate scale free graphs draws back to the "toy" 
model for a system being at a threshold of stability reported in []3T| and to various coherent- 
noise models |52]-[J53] discussed in connection to the standard sandpile model developed 
in self-organized criticality, where the statistics of avalanche sizes and durations also take a 
power law form. The random process in our model begins with no edges at time 0, at a chosen 
vertex i, and adds new edges linked to i, one at a time; each new edge is selected at random, 
uniformly among all outgoing edges not already present. When at some time a performance of 
the system occasionally exceeds a threshold value, the system looses its stability, so that the 
process stops at the given vertex and then starts again at some other vertex lacks of outgoing 
edges. This way the out-degree distribution of the resulting graph equals to the distribution 
of residence times below the threshold which is closely related to the ergodic properties of the 
system and its complexity, |TT| . 

Although it has no a definite relation to the preferential attachment principle |9|], probably, 
both algorithms are similar: we also monitor a global state of the system (while inspecting its 
stability) and introduce an auxiliary function (which is, in some sense, analogous to II (A;)) to 
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characterize the ergodic properties of the system. Nevertheless, the model we propose in the 
present paper has also a distinguishing feature: we allow for the threshold value to quench with 
a given probability rj, so that a natural control parameter in our model is the relative frequency 
of threshold changes. Varying this frequency, one can control the statistics of switching events 
(when the system performance parameter passes the stability threshold) tuning it from the 
truly exponential decay distribution (when the threshold representing the environmental stress 
changes each time a new edge added to the graph) to a power law with some exponent 7 when 
the threshold is frozen. 

Among clear advantages of the proposed algorithm, one can mention its flexibility (the 
possible value of index 7 varies in the wide range (1, 00), the in-degree and out-degree statistics 
can be tuned independently, the variation of frequency rj switches P(k) from an exponential 
decay to a power law) as well as its relative simplicity for both analytical and numerical studies. 

The paper is organized as follows. In Sec. 2, we describe the model we introduce. In Sec. 3, 
we investigate the degree statistics for the random graphs generated in accordance with the 
introduced random procedure. Section 6 is devoted to conclusions. 



2 Motivation and description of the model 

Probably, the most fascinating property of living systems is a surprising degree of tolerance 
against environmental interventions revealing the robustness of the underlying metabolic and 



genetic networks [14, 22, 23]. A typical question which is usually posed in the literature in 



concern to this feature is on the role of the network topology in the error tolerance (see [25 
and references therein). However, in the present work, we bring an opposite question to a 
focus: what is the network topology (say, of the protein-protein interactions map) resulted 
from the very durable evolution process selecting from among species enjoying sufficiently good 
adaptivity? A principle of evolutionary selection of a common large-scale structure of biological 
networks put forward in |23| is closely related to this question. 



The most of biological systems exhibiting a high degree of adaptivity can be referred as 
systems being close to a stability threshold, so that a slight modification of the environment 
would prompt a drastic change in the system performance increasing or reducing the stamina of 
biological spices. From the dynamical systems point of view, such a behavior can be interpreted 
as a forcing of bifurcation parameter through a bifurcation point triggering the system into a 
new state. In our study, we do not refer to any definite physical system displaying such a 



complex behavior. Instead, following pi |, we use a "toy" model, in which the bifurcation 
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parameter (the current state of the system characterized by its net performance) and the 
bifurcation point (the threshold of stability characterizing a challenge of the environments) 
are considered as random independent variables. It is supposed that the system is subjected 
to a bifurcation, when its net performance takes on values above the threshold. It is also 
important that all network's changes acquired in the previous stages of evolution of the graph 
are not discarded, when the system suffers a bifurcation, but still into effect for the forthcoming 
evolution. 

We shall characterize the current performance of the whole system by the real number 
x G [0, 1]. Another real number y G [0, 1] plays the role of the stability threshold. The network 
is supposed to be stable until x < y and is condemned otherwise. We consider i as a random 
variable distributed with respect to some given probability distribution function (pdf) f(u) : 
[0, 1] — > R + . In a deterministic version of the model, / would be an ergodic invariant density 
function on the phase space which has an essential influence on the system dynamics close to 
the stability threshold value. 

The value of threshold y also fluctuates reflecting the changes of environments and is also 
treated as a random variable distributed over the interval [0, 1] with respect to some other pdf 
g(u) : [0, 1] — > R + . In general, / and g are two arbitrary positive integrable functions satisfying 
the normalization condition 



We define the random procedure on the fixed set of N vertices describing the evolution of 
network as a random graph evolving in time. The process begins with no edges at time 0, 
at a chosen vertex %. Given a fixed real number rj G [0, 1], we define a discrete time random 
process in the following way. At time t — 0, the variable x is chosen with respect to pdf /, 
and y is chosen with pdf g. If x < y, we add an edge selected at random, uniformly among all 
edges outgoing from i, the process continues and goes to time t — 1. At time t > 1, one of the 
following events happens: 

i) with probability rj, the random variable x is chosen with pdf / but the threshold y keeps 
the value it had at time t — 1. Otherwise, 

ii) with probability 1 — rj, the random variable x is chosen with pdf /, and the threshold y 
is chosen with pdf g. 

If x > y, the process at the given vertex ends and then resumes at some other vertex chosen 
randomly from the vertices having no outgoing edges yet; if x < y, the process continues at the 
same vertex and goes to time t + 1. 




(2) 
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Eventually, when at some time T limiting the duration of the "monotonous" evolution of 
the graph, the process moves to some other vertex, and i has exactly k = T outgoing edges. 

For example, in Fig. 1 we have presented the results of numerical simulation for the above 
algorithm for a network consisting of 2 x 10 3 vertices, for / = g — 1 and rj — 1. It displays the 
proportion of vertices having degree k versus k in the log-log scale. Circles stay for the out- 
degree k out and diamonds are for the in-degree k in . One can see that the out-degree probability 
distribution is fitted well by a power law, while the in-degree probability distribution forms a 
Gaussian profile centered at the most probable number of incoming links k* n = 11. 

The same algorithm can be extended readily to moderate the incoming edges also. Let us 
consider three random variables x, y, and z that are the real numbers distributed in accordance 
to the distributions /, g, and v within the unit interval [0, 1]. We assume that x represents the 
current performance of the system, y and z are the thresholds for emission and absorption of 
edges respectively. 

The generalized process begins on the set of N vertices with no edges at time 0, at a chosen 
vertex %. Given two fixed numbers i] £ [0, 1] and v £ [0, 1], the variable x is chosen with respect 
to pdf /, y is chosen with pdf g, and z is chosen with pdf v, we draw edge outgoing from 
% vertex and entering j vertex if x < y and x < z and continue the process to time t — 1. 
Otherwise, if x > y (x > z), the process moves to other vertices having no outgoing (incoming) 
links yet. 

At time t > 1, one of the three events happens: 

i) with probability 77, the random variable x is chosen with pdf / but the thresholds y and 
z keep their values they had at time t — 1. 

ii) with probability 1 — 77, the random variable x is chosen with pdf /, and the thresholds y 
and z are chosen with pdf g and v respectively. 

iii) with probability is, the random variable x is chosen with pdf /, and the threshold z is 
chosen with pdf v but the threshold y keeps the value it had at time t — 1. 

If x > y, the process stops at % vertex and then starts at some other vertex having no 
outgoing edges yet. If x > z, the accepting vertex j is blocked and does not admit any more 
incoming link (provided it has any). If x < y and x < z, the process continues at the same 
vertex % and goes to time t + 1. 

Figure 2 displays an example of a graph generated in accordance this algorithm on the set 
of 100 vertices for the case of uniform densities f — g — v — 1 for rj = v = 1. In Fig. 3, we have 
presented the distribution P(k) vs. k in the log-log scale over 10 5 vertices. Here, the circles 
stay for outgoing degrees, and diamonds are for incoming degrees. For k 3> 1, both profiles 
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enjoy a power law decay with 7 in = 7 out = 2. 

In the forthcoming section, we demonstrate that for some probability densities /, g, and v 
the resulting probability distribution P(k) meets asymptotically a power law (|l]) as k 3> 1 
where the index 7 can vary in the wide range 7 G (l,oo). Moreover, by means of adjusting 
probability densities f,g, and v from a proper class of functions, one can tune 7i n and 7 out to 
different values. 

If v = 0, both thresholds y and z have synchronized dynamics, and sliding the value of rj 
form to 1, one can tune the statistics of out-degrees and in-degrees simultaneously out from 
the pure exponential decay to the power laws provided /, g, and v belong to the proper class 
of functions. When v — 1, the threshold z changes each time but the threshold y is frozen. In 
this case, P(h a ) decays exponentially for any choice of v, but P(k out ) still meets a power law 
for k out > 1. 



3 Probability degree distribution P(k) 



Below, we compute the probability P(k) that a given vertex has precisely k outgoing edges 
as a result of the random procedure defined in the previous section. For sake of clarity, we 
investigate the statistics of switching events quenching the monotonous evolution of the graph 
in the simplest case when only one threshold (y) is present. All computations can be extended 
readily to the general case when both thresholds are present. 

Directly from the definitions of Sec. 2, one can show that 



dy g{y) / dx f(x) 



P(0) 



In t > 1 time steps, the system can either "survive" (S) or "die" (D). Both issues can take 
place either in the "correlated" way (with probability r), see i) (we denote these scenarios as 
S c and D c consequently), or in the " uncorrelated" way (with probability 1 — 77, see ii) (S u and 
D u , consequently). For example, for k = 1, we have 

P(l) = P[SD C ] + P[SD U ] = V B{1) + (1 - r))A 2 {l) 



where we have defined, for n = 0, 1, 2, ... , 



A{n) 



dy g{y) 



dx f(x) 



U0 



and 



Bin) 



dy g{y) 



dx f{x) 



Uo 



du f(u) = A{n) - A(n + 1). 



(3) 



(4) 
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Similarly, 

P(2) = P[SS C D L ] + P[SS C D U ] + P[SS U D C ] + P[SS U D U ], 

where 

P[SS C D C ] = V 2 B(2), P[SS C D U ] = 77(1 - rj)A(2)B(0), 

P[SS U D C ] = 77(1 - 7j)A(l)fl(l), PfS^Aj = (1 - ?7) 2 A 2 (1)S(0). 

The general formula for k > 3 is tedious and can be found in ||31|| . 

Computations can be essentially simplified with the use of generating functions. Let us 
define the generating function for P(k) by 

1 d k P(s 



k=0 



(5) 



s=0 



and introduce the new auxiliary functions, 

p(l)=r} l A{l + i), for/>l, p(0) = 0, 
g(/) = (l-r / )^- 1 (l), for/>l, g(0) = 0, 
r (Z) =^[r/ J B(/ + 1) + (1-t/)A(Z + 1)5(0)], for / > 1, r(0) = 0, 

p = 77^(1) + (1-77)A(1)S(0). (6) 
Then, the use of convolutions property []35| yields 

P(s) = B(0) + ps + [r(s) + pp(s)q(s) + P A(l)q(s) + A(l)q(s)f(s)} , (7) 

where p(s), q(s), r(s) are the generating functions of p(l), q(l),r(l), respectively. 

In the marginal cases of rj = and rj = 1, the probability P(k) can be readily calculated. 

3.1 The uncorrelated case, 77 = 

This case corresponds to the external stress being sensitive to any modification of a network 
and changing its figure as soon as a new edge appears in the graph. For 77 = 0, one has 
p(s) = f(s) = 0, and q(s) = s(l — sA(l)) -1 , p = A(1)B(0), then equation ([?D gives 

Applying the inverse formula in (|5|) to equation (^j, one arrives at 



P v=o (k)=A k (l)B(0) 



1 rv 

dy g{y) \ dx f(x) 

JO 



k r i „i 

dy g{y) \ dx f(x) 

Jy 



Therefore, in the "uncorrelated" case (77 = 0), for any choice of the pdf f{u) and g(u), the 
probability P(k) decays exponentially. 
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3.2 The correlated case, rj — 1 



Another marginal case corresponds to the external stress remaining unchanged during all 
time of network evolution. For rj = 1, one has q(s) = 0, p = B(l), and 



sr 



s) = s n+1 B(n + 1) = J2 sUB[ 



n) 



n=l 



n=2 



so that formula (J?]) yields 
Consequently, 



P v=1 (s) = B(s) 



P n= i(k)=B(k) 



dy g{y) 



dx f(x) 



du f(u). 



(9) 



(10) 



11' 



In the "correlated" case, r\ — 1, many different types of behaviour are possible, depending on 
the form of the pdf f(u) and g(u). We will examine an important class of f(u) and g(u), for 
which P v= i(k) can be explicitly computed from equation fllpp . We will take 

f(u) = (1 + a)u a , a > -1, 

g(u) = (l + P)(l-uy, 0>-l. 

Equation fllCf ) gives in this case: 

_ r(2 + g r(i + fc(i + a)) r(2 + g r(i + (fc + i)(i + a)) 

77-11 J T(2 + /5 + fc(l + a)) r(2 + p + (k + 1)(1 + a)) ' 

Using Stirling's approximation, we get for large k » 1: 



;i + /3) r(2 + /?) (i + q 



,-1-/3 



For different values of (3, the exponent of the threshold distribution, we get all possible power 
law decays of P v =i{k). Notice that the exponent (—2 — 0) characterizing the decay of P v= i(k) 
is independent of the distribution f(u) of the state variable. 

Therefore, for rj = 1, one has to provide just one power law function (which defines the 
statistics of the threshold) to achieve a power law for P(k) characterized by the exponent 

7 e (i,oo). 

3.3 The case of uniform densities / = g = 1 

For the case of uniform densities / = g = 1, and after some tedious but trivial computation, 
we get from equation (|7j): 



P(s) 



l + (l-77) M (s) 



(12) 
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where /x(s) is defined by 

ln(l - 17s) 

fi(s) = . 

rjs 

The asymptotic behavior of P(k) for k large is determined by the singularity of the generating 
function P(s) that is closest to the origin. 

For 77 = 0, the generating function has a simple pole P(s) = (2 — s) , and therefore P(k) 
decays exponentially, which agrees with the general conclusion of Sec. 3.1. 

For intermediate values 1 > 77 > 0, the generating function P(s) has two singularities. One 
pole, s = Sq, corresponds to the vanishing denominator 1 + (1 — 77)7/(5), where So is the unique 
nontrivial solution of the equation 

- ln(l - 77s) = s—^— . (13) 

1 — 7] 

Another singularity, s = S\, corresponds to the vanishing argument of the logarithm, s% = f]^ 1 , 
such that Si > s > 1. The dominant singularity of P(s) is of polar type, and the corresponding 
decay of P(k) is asymptotically exponential for the degrees k larger than k c ~ ln(so) -1 , with a 
rate ln(so) that vanishes like 1 — r] as 7] tends to 1. 

When 7] tends to one, two singularities sq and s± merge. More precisely, we have 

P v =i{s) = 2 . (14) 

The corresponding dominant term in (|T^ ) is of order 0(&;~ 2 ) ^j. This obviously agrees with 
the exact result one can get from equation fllCf ), with f(x) — g(x) = 1: 

P ^ [k) = (fc +!)'(* + 2) ' <15) 
In the case of uniform densities, it is possible to get an expression of P(k) for all k, and for any 
value of 77, by applying the inversion formula @ to equation ([PD and use of the formulae [j36 
Combining all this, we get 

T] x ^ 7] / 1 — 77 



where c(m, Z) is defined by 
c(m, I) = ml 



/l+«2 + --- + «m =« 
ii>l 



(h + h) (h + 1) (I ~ h ~ h) ■ ■ ■ (lm-1 + l)(l-h-l 2 l m -l) (lm + 1) 
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When rj 7^ 0, there is an alternative way of writing the previous expression: 



m - <rrfe + E (t ._ i+ V- i+2) X> - , (17) 

where 6(m, /) is defined by 
b(m,l) = 



(ii + 1) (z 2 + 1) ■ • • (i m _i + 1) (z m + 1) ' 

h+i2-\ Hi =1 



In Fig. 4, we have plotted these asymptotes in the log-linear scale for the consequent values 
77 = 0, 77 = 0.5, 77 = 0.7, 77 = 0.9, T] = 1 (bottom to top) together with the curves representing 
analytical results ([T5| - p~6| ) . 

Notice that for rj ^ 1, we only plot distributions -P(fc) up to relatively small degrees (k = 16, 
A; = 25) correspondent to short quiescent times in the system, since these times are already 
bigger then the crossover value kjrj) ~ l/ln(s ) to the exponential decay exp (— ln(s ) k) 
(where sq is defined by the Eq. fll3[)). For much longer times, very few survivals are observed, 
and the statistics gets bad. Of course, k c (r]) grows as the parameter 77 tends to 1, so that we 
have good statistics for larger k. 



4 Discussion and Conclusion 

In this paper, we have presented a simple random procedure inducing various types of 
random graphs and the scale free random graphs among others. The model we use essentially 
differs from the preferential attachment approach || discussed in the literature before. The 
proposed approach is simple and more flexible both for numerical simulations and analytical 
studies. The model draws back to the original "toy" model for a system being at a threshold 
of stability 

In the proposed model, the degree statistics of a random graph depends from the value of 
the control parameter 77 stirring the pure exponential statistics of degrees (at 77 = 0, when the 
threshold value is sensitive to each new edge risen in the graph) to a power law (at 77 = 1, when 
the threshold value is frozen). The exponent characterizing the power law depends as 7 = 2 + (3 
from the power law probability density distribution of the threshold g(u) oc (1 — u)^, (3 > — 1, 
and can vary in the wide range 7 G (1, 00). For the intermediate values of 77, the decay rate is 
mixed. In particular, we showed that, even if the asymptotic behavior is exponential, when r\ 
is close to 1 this only holds for very large degrees k of order of (1 — rf)^ 1 . The model can be 
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readily generalized to various particular types of networks. The crucial advantage of the model 
is that, for most physically meaningful distributions, it can be solved analytically. In other 
numerical solution can be readily found. 
In the present paper, we have limited the consideration to the special cases of uniform and 
power law distributions of thresholds. Indeed, one can take other (normalized) distributions 
for /, g, and v in this model and obtain different decaying asymptotic profiles for P(k) either 
for synchronous v = or asynchronous v — 1 dynamics of thresholds, they have, in general, 
nothing to do with power laws at r] — 1, but still uniformly exponential at r] = 0. 
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Figure 1: The proportion of vertices of degree k versus k in the log-log scale (N = 2000, / = g = 1, 77 = 1). 
Out-degree statistics (circles) satisfies a power law, while the in-degree probability distribution is Gaussian with 
the most probable number of incoming links k* n = 11. 
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Figure 3: The distribution P(k) vs. k in the log-log scale generated on 10 5 vertices, / = g = v = 1, 77 = 1. 
Here, the circles stay for outgoing degrees and diamonds are for incoming degrees. Both profiles enjoy a power 
law decay with 7 = 2. 
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Figure 4: The distribution P(k) in log-linear scale for different values of the parameter 77, rj = 0., r\ = 0.5, 
rj = 0.7, rj = 0.9 and 77 = 1.0 (bottom to top). Straight lines corresponds to the pure exponential decay 
(~ 2~ fe ) observed in the case 77 = 0, the results obtained from ( |l6| ) for the intermediate values of rj, and to the 
result (Eq) for 77 = 1.0 (bottom to top). 



